Load libraries

library(data.table)
library(tictoc)
library(ggplot2)
library(plyr)
library(dplyr)
library(signal)
library(roll)
library(pracma)
library(plotly)
source("peakdet.R")
library(rootSolve)
library(scales)
library(tidyr)
library(qpcR)

library(DT)

Import data (shiny version uses input from ui)

BoatdfAl <- data.frame(fread(choose.files(default = "C:/Users/aaron.beach/CaptureU", caption = "Select Boat File")))
PelvisdfAl <- data.frame(fread(choose.files(default = "C:/Users/aaron.beach/CaptureU", caption = "Select Pelvis File")))
ThoraxdfAl <- data.frame(fread(choose.files(default = "C:/Users/aaron.beach/CaptureU", caption = "Select Thorax File")))
MinimaxdfAl <- data.frame(fread(choose.files(default = "C:/Users/aaron.beach/CaptureU", caption = "Select Minimax File")))

Calculate Sample Rate, then trim data to between 60-90 seconds

t = length(BoatdfAl$time_s)-1
SampleRateTS <- matrix(0, ncol = 1, nrow = t)
for (i in 1:t) {
  SampleRateTS[i,] <- 1/(BoatdfAl$time_s[i+1]-BoatdfAl$time_s[i])
}
SampleRate <- median(SampleRateTS)



start = SampleRate*60
end = SampleRate*90

BoatdfAl <- BoatdfAl[start:end,]
PelvisdfAl <- PelvisdfAl[start:end,]
ThoraxdfAl <- ThoraxdfAl[start:end,]

Reset the index/column V1 from zero

row.names(BoatdfAl)<-1:nrow(BoatdfAl)
row.names(PelvisdfAl)<-1:nrow(PelvisdfAl)
row.names(ThoraxdfAl)<-1:nrow(ThoraxdfAl)

BoatdfAl$V1<-1:nrow(BoatdfAl)
PelvisdfAl$V1<-1:nrow(PelvisdfAl)
ThoraxdfAl$V1<-1:nrow(ThoraxdfAl)

Adjust values i.e. if using the aligned files, the gyro is converted to rad/s, so to convert back to deg/s multiply by 57.2958

The global angle data is adjusted based on the offset from the boat orientation during the static pose. Done manually but developing an auto way to do this.

#PelvisdfAl$gx_deg.s = PelvisdfAl$gx_deg.s*57.2958
#ThoraxdfAl$gx_deg.s = ThoraxdfAl$gx_deg.s*57.2958
PelvisdfAl$z_deg2 = PelvisdfAl$z_deg-BoatdfAl$z_deg+10
ThoraxdfAl$z_deg2 = ThoraxdfAl$z_deg-BoatdfAl$z_deg+8

Create spectrum plots to identify filtering frequency. Filter each variable using a butterworth filter.

spectrum(BoatdfAl$"ax_m.s.s", method = "ar")

bf <- butter(4, 0.1, type="low")
BoatdfAl$AccFilt <- filtfilt(bf, BoatdfAl$"ax_m.s.s")


#spectrum(PelvisdfAl$"gx_deg.s", method = "ar")

bf <- butter(4, 0.1, type="low")
PelvisdfAl$gyroxfilt <- filtfilt(bf, PelvisdfAl$"gx_deg.s")

#spectrum(ThoraxdfAl$"gx_deg.s", method = "ar")

bf <- butter(4, 0.1, type="low")
ThoraxdfAl$gyroxfilt <- filtfilt(bf, ThoraxdfAl$"gx_deg.s")

#spectrum(PelvisdfAl$z_deg2, method = "ar")

bf <- butter(4, 0.1, type="low")
PelvisdfAl$z_degfilt <- filtfilt(bf, PelvisdfAl$z_deg2)

#spectrum(ThoraxdfAl$z_deg2, method = "ar")

bf <- butter(4, 0.1, type="low")
ThoraxdfAl$z_degfilt <- filtfilt(bf, ThoraxdfAl$z_deg2)
#spectrum(ThoraxdfAl$z_degfilt, method = "ar")

To split the data into strokes:

Identify acceleration peaks (Event) Use the second column (first and 3rd column are the values for where the “mountain” starts and ends) Check the plot to see if these were identified correctly, if not adjust findpeaks criteria

Event <- findpeaks((BoatdfAl$AccFilt), minpeakheight = 1, minpeakdistance = 50)

ev <-  ggplot() + geom_line(aes(x=BoatdfAl$V1, y=(BoatdfAl$AccFilt))) + geom_point(aes(x=Event[,2], y=Event[,1]))
ggplotly(ev)


Event <- Event[,2]
Event <- sort(Event, decreasing = FALSE)

If the first peak doesn’t have a zero value prior (due to where the data was trimmed), remove it.

if(which.min(BoatdfAl$AccFilt[1:Event[1]]) >0){
  Event <- Event[-1]
}

Create stroke cycle start/end as the zero prior to a peak.

zero <- NULL
cycle <- NULL
cycle[1] = tail((which(BoatdfAl$AccFilt[1:Event[1]] <0)), n=1)


for (r in 1:(length(Event))){
  t = r+1
  # framer = Event[r,2] - BoatdfAl$V1[1]
  # framet = Event[t,2] - BoatdfAl$V1[1]
  
  tryCatch((cycle[t] = tail((which(BoatdfAl$AccFilt[Event[r]:Event[t]] <0)), n=1) + Event[r]),error = function(e) print(NA))
}

Cut up data into stroke cycles, i.e. cycle 1:2, cycle 2:3 etc.

This code often throws up errors due to uneven rows, NAs etc. Not an issue in RStudio, but in RShiny the app stops. tryCatch ignores the error.

BoatAcc = NULL
BoatAcc1 = NULL
BoatAcc2 = NULL
PelvisGyro = NULL
PelvisGyro1 = NULL
PelvisGyro2 = NULL
ThoraxGyro = NULL
ThoraxGyro1 = NULL
ThoraxGyro2 = NULL
PelvisRot = NULL
PelvisRot1 = NULL
PelvisRot2 = NULL
ThoraxRot = NULL
ThoraxRot1 = NULL
ThoraxRot2 = NULL

for (c in 1:(length(Event))){
  t = c+1
  tryCatch((BoatAcc[[c]] = BoatdfAl$AccFilt[cycle[c]:cycle[t]]), error = function(e) print(NA))
  tryCatch((PelvisGyro[[c]] = PelvisdfAl$gyroxfilt[cycle[c]:cycle[t]]),error = function(e) print(NA))
  tryCatch((ThoraxGyro[[c]] = ThoraxdfAl$gyroxfilt[cycle[c]:cycle[t]]),error = function(e) print(NA))
  tryCatch((PelvisRot[[c]] = PelvisdfAl$z_degfilt[cycle[c]:cycle[t]]),error = function(e) print(NA))
  tryCatch((ThoraxRot[[c]] = ThoraxdfAl$z_degfilt[cycle[c]:cycle[t]]),error = function(e) print(NA))
  
}

Use Thorax Gyro to determine if Left or Right is the first stroke. If first stroke max is positive, then it is a Left Stroke.

if (max(ThoraxRot[[1]]) >0){
  Side1 <- "Left"
  
}else{
  Side1 <- "Right"
  
}

Cut up variables per stroke. Label Left and Right according to the above.

for (i in seq(2,length(BoatAcc),2)){
  k = i-1
  tryCatch((BoatAcc2[[i-1]] = BoatAcc[[i]]), error = function(e) print(NA))
  tryCatch((BoatAcc1[[i-1]] = BoatAcc[[k]]), error = function(e) print(NA))
}

BoatAcc1 <- plyr::compact(BoatAcc1)
BoatAcc2 <- plyr::compact(BoatAcc2)

if (Side1 == "Left"){
  BoatAccLeft <- BoatAcc1
    BoatAccRight <- BoatAcc2
} else {
    BoatAccRight <- BoatAcc1
    BoatAccLeft <- BoatAcc2
}



for (i in seq(2,length(PelvisGyro),2)){
  k = i-1
  tryCatch((PelvisGyro2[[i-1]] = PelvisGyro[[i]]), error = function(e) print(NA))
  tryCatch((PelvisGyro1[[i-1]] = PelvisGyro[[k]]), error = function(e) print(NA))
}

PelvisGyro1 <- plyr::compact(PelvisGyro1)
PelvisGyro2 <- plyr::compact(PelvisGyro2)

if (Side1 == "Left"){
  PelvisGyroLeft <- PelvisGyro1
    PelvisGyroRight <- PelvisGyro2
} else {
    PelvisGyroRight <- PelvisGyro1
    PelvisGyroLeft <- PelvisGyro2
}

for (i in seq(2,length(ThoraxGyro),2)){
  k = i-1
  tryCatch((ThoraxGyro2[[i-1]] = ThoraxGyro[[i]]), error = function(e) print(NA))
  tryCatch((ThoraxGyro1[[i-1]] = ThoraxGyro[[k]]), error = function(e) print(NA))
}


ThoraxGyro1 <- plyr::compact(ThoraxGyro1)
ThoraxGyro2 <- plyr::compact(ThoraxGyro2)

if (Side1 == "Left"){
  ThoraxGyroLeft <- ThoraxGyro1
    ThoraxGyroRight <- ThoraxGyro2
} else {
    ThoraxGyroRight <- ThoraxGyro1
    ThoraxGyroLeft <- ThoraxGyro2
}

for (i in seq(2,length(PelvisRot),2)){
  k = i-1
  tryCatch((PelvisRot2[[i-1]] = PelvisRot[[i]]), error = function(e) print(NA))
  tryCatch((PelvisRot1[[i-1]] = PelvisRot[[k]]), error = function(e) print(NA))
}

PelvisRot1 <- plyr::compact(PelvisRot1)
PelvisRot2 <- plyr::compact(PelvisRot2)

if (Side1 == "Left"){
  PelvisRotLeft <- PelvisRot1
    PelvisRotRight <- PelvisRot2
} else {
    PelvisRotRight <- PelvisRot1
    PelvisRotLeft <- PelvisRot2
}

for (i in seq(2,length(ThoraxRot),2)){
  k = i-1
  tryCatch((ThoraxRot2[[i-1]] = ThoraxRot[[i]]), error = function(e) print(NA))
  tryCatch((ThoraxRot1[[i-1]] = ThoraxRot[[k]]), error = function(e) print(NA))
}

ThoraxRot1 <- plyr::compact(ThoraxRot1)
ThoraxRot2 <- plyr::compact(ThoraxRot2)

if (Side1 == "Left"){
  ThoraxRotLeft <- ThoraxRot1
    ThoraxRotRight <- ThoraxRot2
} else {
    ThoraxRotRight <- ThoraxRot1
    ThoraxRotLeft <- ThoraxRot2
}

Minimax data is 100Hz and not synced. So deal with it separately. Divide into the same period of time, using time column (rather than row index).

Minimaxdf2 <-  subset(MinimaxdfAl, `Time`>=60 & `Time`<=90)

Minimax has a column “stks” which identifies the start of each stroke, alternating 1 and 2, with blanks in between (imported as NA).

stks1 = which(Minimaxdf2$Stk == 1)
stks2 = which(Minimaxdf2$Stk == 2)
stks = c(stks1, stks2)
stks = sort(stks)

Cut up the Minimax Velocity column per Stroke.

MinimaxdfVelocity = NULL
for (i in 2:length(stks)){
  k = i-1
  
  MinimaxdfVelocity[[k]] = Minimaxdf2[stks[k]:stks[i],"VAcc"]
  
}

MinimaxdfVelocity1=NULL
MinimaxdfVelocity2=NULL

for (i in seq(2,length(MinimaxdfVelocity),2)){
  k = i-1
  tryCatch((MinimaxdfVelocity2[[i-1]] = MinimaxdfVelocity[[i]]), error = function(e) print(NA))
  tryCatch((MinimaxdfVelocity1[[i-1]] = MinimaxdfVelocity[[k]]), error = function(e) print(NA))
  }

MinimaxdfVelocity1 <- plyr::compact(MinimaxdfVelocity1)
MinimaxdfVelocity2 <- plyr::compact(MinimaxdfVelocity2)

In a similar process to above, whether the first stroke is a left or right is determined by the Minimax Roll column. If the peak roll for the first stroke is positive, then the first stroke is Right.

MinimaxdfRoll = NULL
for (i in 2:length(stks)){
  k = i-1
  MinimaxdfRoll[[k]] = Minimaxdf2[stks[k]:stks[i],"Roll"]
  
}


if (max(MinimaxdfRoll[[1]]) >2){
  side = "right"
  MinimaxdfVelocityRight <- MinimaxdfVelocity1
  MinimaxdfVelocityLeft <- MinimaxdfVelocity2
  
}else{
  side = "left"
  MinimaxdfVelocityLeft <- MinimaxdfVelocity1
  MinimaxdfVelocityRight <- MinimaxdfVelocity2
  
}

MinimaxdfVelocityRight <- plyr::compact(MinimaxdfVelocity1)
MinimaxdfVelocityLeft <- plyr::compact(MinimaxdfVelocity2)

This process is common for all variables, normalising each stroke with varying data points, into a 0-100% scale.

MinimaxdfVelocityRightNorm = NULL
for (i in 1:length(MinimaxdfVelocityRight)){
  xvalues = 1:length(MinimaxdfVelocityRight[[i]])
  yvalues = MinimaxdfVelocityRight[[i]]
  output_x_vals <- seq(1,length(xvalues),length.out=101)
  #
  # compute the interpolated values; this would be done for each input time series
  #
  #interp_output[[,i]]<- 
  
  MinimaxdfVelocityRightNorm[[i]] <- approx(x=xvalues, y=yvalues, xout=output_x_vals)$y
}

MinimaxdfVelocityLeftNorm = NULL
for (i in 1:length(MinimaxdfVelocityLeft)){
  xvalues = 1:length(MinimaxdfVelocityLeft[[i]])
  yvalues = MinimaxdfVelocityLeft[[i]]
  output_x_vals <- seq(1,length(xvalues),length.out=101)
  #
  # compute the interpolated values; this would be done for each input time series
  #
  #interp_output[[,i]]<- 
  
  MinimaxdfVelocityLeftNorm[[i]] <- approx(x=xvalues, y=yvalues, xout=output_x_vals)$y
}

Mean, Standard deviation, positive and negative error bars, and max (magnitude and %Stroke) are also calculated.

MinimaxdfVelocityRightNormCombined = do.call(cbind, MinimaxdfVelocityRightNorm)
MinimaxdfVelocityRightNormMean = rowMeans(MinimaxdfVelocityRightNormCombined)
MinimaxdfVelocityRightNormSD = apply(MinimaxdfVelocityRightNormCombined[,-1], 1, sd)
MinimaxdfVelocityRightNormPEB = MinimaxdfVelocityRightNormMean + MinimaxdfVelocityRightNormSD
MinimaxdfVelocityRightNormNEB = MinimaxdfVelocityRightNormMean - MinimaxdfVelocityRightNormSD

MinimaxdfVelocityRightNormCombined2 = as.data.frame(MinimaxdfVelocityRightNormCombined)
MinimaxdfVelocityRightNormCombined2$perc = c(0:100)
MinimaxdfVelocityRightNormCombined2 <- melt(MinimaxdfVelocityRightNormCombined2, id.vars="perc")
## Warning in melt(MinimaxdfVelocityRightNormCombined2, id.vars = "perc"): The melt
## generic in data.table has been passed a data.frame and will attempt to redirect
## to the relevant reshape2 method; please note that reshape2 is deprecated, and
## this redirection is now deprecated as well. To continue using melt methods from
## reshape2 while both libraries are attached, e.g. melt.list, you can prepend the
## namespace like reshape2::melt(MinimaxdfVelocityRightNormCombined2). In the next
## version, this warning will become an error.
MinimaxdfVelocityRightNormMeanMax <- round(max(MinimaxdfVelocityRightNormMean),2)
MinimaxdfVelocityRightNormMeanMaxPC <- which.max(MinimaxdfVelocityRightNormMean)-1

All strokes are plotted together.

#output$MinimaxVelRight<- renderPlotly({
  
  fig <- plot_ly(MinimaxdfVelocityRightNormCombined2, x = ~perc, y = ~value, showlegend = FALSE) %>% add_lines(color = ~ordered(variable))  %>% 
    layout(xaxis = list(title = "% of Stroke"), yaxis = list(title = "Minimax Velocity (m/s)"))
  fig
## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
#})

Process is repeated on the left side

Both left and right mean and error bars are plotted together.

#output$MinimaxVelCombined<- renderPlotly({
  
  fig <- plot_ly(x = c(0:100)) %>%
    add_trace(y = MinimaxdfVelocityRightNormMean, type = 'scatter', mode = 'lines',
              line = list(color = 'rgba(250, 0, 0, 1)'),
              showlegend = TRUE, name = 'Right Stroke Minimax Velocity')%>% 
    add_trace(y = MinimaxdfVelocityRightNormPEB, type = 'scatter', mode = 'lines',
              fill = 'tonexty', fillcolor='rgba(255, 10, 10, 0.2)', line = list(color = 'rgba(0,0,0,0)'),
              showlegend = FALSE, name = 'Lowerlimit')%>% 
    add_trace(y = MinimaxdfVelocityRightNormNEB, type = 'scatter', mode = 'lines',
              fill = 'tonexty', fillcolor='rgba(255, 10, 10, 0.2)', line = list(color = 'rgba(0,0,0,0)'),
              showlegend = FALSE, name = 'Lowerlimit')%>% 
    add_trace(y = MinimaxdfVelocityLeftNormMean, type = 'scatter', mode = 'lines', line = list(color = 'green'),
              showlegend = TRUE, name = 'Left Stroke Minimax Velocity')%>% 
    add_trace(y = MinimaxdfVelocityLeftNormPEB, type = 'scatter', mode = 'lines',
              fill = 'tonexty', fillcolor='rgba(0,100,80,0.2)', line = list(color = 'rgba(0,0,0,0)'),
              showlegend = FALSE, name = 'Lowerlimit')%>% 
    add_trace(y = MinimaxdfVelocityLeftNormNEB, type = 'scatter', mode = 'lines',
              fill = 'tonexty', fillcolor='rgba(0,100,80,0.2)', line = list(color = 'rgba(0,0,0,0)'),
              showlegend = FALSE, name = 'Lowerlimit')%>% 
    layout(         xaxis = list(title = "% of Stroke",
                                 showgrid = TRUE,
                                 showline = FALSE,
                                 showticklabels = TRUE,
                                 ticks = 'outside',
                                 zeroline = TRUE),
                    yaxis = list(title = "Minimax Velocity (m/s)",
                                 showgrid = TRUE,
                                 showline = FALSE,
                                 showticklabels = TRUE,
                                 ticks = 'outside',
                                 zeroline = TRUE),
                    legend = list(x = 0.5, y = 0.9)) %>%
    config(displayModeBar = F)
  
  
  fig
  #})

The same process is repeated for the Boat Acceleration.

  fig
#})

Pelvis Gyro (Rotational Velocity)

  fig
#})

Thorax Gyro (Rotational Velocity)

  fig
#})

All means are plotted together for the right side

  fig
#})

And Left side.

  fig
#})

Max values are collated into a datatable.

SequencedataRight <- data.frame(MinimaxdfVelocityRightNormMeanMax, MinimaxdfVelocityRightNormMeanMaxPC, PelvisGyroRightNormMeanMax,  PelvisGyroRightNormMeanMaxPC, ThoraxGyroRightNormMeanMax, ThoraxGyroRightNormMeanMaxPC)
colnames(SequencedataRight) <- c("Minimax Velocity Max (m/s)", "%", "Pelvis Rotation Velocity Max (deg/s)",  "%", "Thorax Rotation Velocity Max (deg/s)", "%")
SequencedataLeft <- data.frame(MinimaxdfVelocityLeftNormMeanMax, MinimaxdfVelocityLeftNormMeanMaxPC, PelvisGyroLeftNormMeanMax,  PelvisGyroLeftNormMeanMaxPC, ThoraxGyroLeftNormMeanMax, ThoraxGyroLeftNormMeanMaxPC)
colnames(SequencedataLeft) <- c("Minimax Velocity Max (m/s)", "%", "Pelvis Rotation Velocity Max (deg/s)",  "%", "Thorax Rotation Velocity Max (deg/s)", "%")

Sequencedata <- rbind(SequencedataRight,SequencedataLeft)


#output$table_VelSequence <- renderDataTable({ 
  
  
  datatable(Sequencedata,  rownames = c("Right", "Left"), options = list(columnDefs = list(list(className = 'dt-center', targets = "_all")),dom='t', ordering=F))%>% formatStyle(
    columns = c(3,4),
    #c('250m splits (secs)', "Race 2", "Avg Velocity (m/s)", "Avg Prog Speed (%)"),
    valueColumns = 0,
    target = 'cell',
    backgroundColor = 'rgba(112, 165, 255, 1)')
#})

Repeat for Pelvis orientation (global angle)

  fig
#})

And Thorax orientation (global angle)

  fig
#})

Plot all the means together for the right side.

  fig
#})

And the left side.

  fig
#})

And collate the means into a datatable.

  datatable(SequenceAngledata,  rownames = c("Right", "Left"), options = list(columnDefs = list(list(className = 'dt-center', targets = "_all")),dom='t', ordering=F))%>% formatStyle(
    columns = c(3,4),
    #c('250m splits (secs)', "Race 2", "Avg Velocity (m/s)", "Avg Prog Speed (%)"),
    valueColumns = 0,
    target = 'cell',
    backgroundColor = 'rgba(112, 165, 255, 1)')
#})